#setwd("filepath")
path<-getwd()
library(coda)

load("ecrt_demsq_models.rda")
load("ecrt_dif_models.rda")
load("ecrt_data_NAs.rda")

ecrt.models<-the.hi.model.list
ecrt.dif.models<-the.hi.model.dif.list

ecrt.dat<-newest.dat
dat<-na.omit(newest.dat)
dat$end<-ifelse(dat$year==1998,1,0)
dat$end<-ifelse(dat$ecrt==1,1,dat$end)
dat$end[dat$ccode==265&dat$year==1989]<-1
dat$id<-rev(cumsum(rev(dat$end)))
dat<-dat[order(dat$id, dat$year),]

load("e_ccpr_demsq_models.rda")
load("e_ccpr_dif_models.rda")
load("e_ccpr_data_NAs.rda")

ccpr.models<-the.hi.model.list
ccpr.dif.models<-the.hi.model.dif.list

### extract coefficients for E Court models

ests.list<-list()
for (i in 1:10){
	ests1<-as.matrix(ecrt.models[[i]][[1]])	
	ests2<-as.matrix(ecrt.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.list[[i]]<-ests
	}

ests.dif.list<-list()
for (i in 1:50){
	ests1<-as.matrix(ecrt.dif.models[[i]][[1]])	
	ests2<-as.matrix(ecrt.dif.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.dif.list[[i]]<-ests
	}

#### effect of democracy on E Court from democracy squared models
probs.list<-list()
dem<-seq(from=-1.56, to=2.11, by=.01)
for (i in 1:10){

	probs<-matrix(nrow=368, ncol=4)

	for (j in 1:368){
		pr<-
		plogis(
		ests.list[[i]][,46]+
		ests.list[[i]][,42]*dem[j]+
		ests.list[[i]][,44]*dem[j]^2+
		ests.list[[i]][,45]*20
		)
	
		probs[j,1]<-quantile(pr, probs=.05)
		probs[j,2]<-mean(pr)
		probs[j,3]<-quantile(pr, probs=.95)
		probs[j,4]<-dem[i]
		}

	probs.list[[i]]<-probs

	}

### extract estimates for CCPR models

ests.list.2<-list()
for (i in 1:10){
	ests1<-as.matrix(ccpr.models[[i]][[1]])	
	ests2<-as.matrix(ccpr.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.list.2[[i]]<-ests
	}

ests.dif.list.2<-list()
for (i in 1:50){
	ests1<-as.matrix(ccpr.dif.models[[i]][[1]])	
	ests2<-as.matrix(ccpr.dif.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.dif.list.2[[i]]<-ests
	}

### effect of democracy on CCPR from democracy squared models
probs.list.2<-list()
dem<-seq(from=-1.56, to=2.11, by=.01)
for (i in 1:10){

	probs<-matrix(nrow=368, ncol=4)

	for (j in 1:368){
		pr<-
		plogis(
		ests.list.2[[i]][,35]+
		ests.list.2[[i]][,31]*dem[j]+
		ests.list.2[[i]][,33]*dem[j]^2+
		ests.list.2[[i]][,34]*20
		)
	
		probs[j,1]<-quantile(pr, probs=.05)
		probs[j,2]<-mean(pr)
		probs[j,3]<-quantile(pr, probs=.95)
		probs[j,4]<-dem[i]
		}

	probs.list.2[[i]]<-probs

	}

tr.gr<-rgb(150,150,150,150,maxColorValue=255)

pdf("ecrt_demsq.pdf", width=11, height=5)

#dev.new(width=11,height=5)
par(mfrow=c(1,2), mar=c(4, 4, 2, 1), cex.lab=1.25,cex.main=1.5)

plot(density(ecrt.dat$mean), ylim=c(0, 1), xlim=c(-2, 2.5), main="European Court", xlab="Latent Democracy", ylab="Pr(Acceptance)", axes=F, type="n")
abline(h=c(0,0.2,0.4,0.6,0.8,1), col=tr.gr, lty=3, lwd=1.5)
polygon(density(ecrt.dat$mean), col="grey")
axis(1, at=c(-2, -1, 0, 1, 2))
axis(2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1))
polygon(density(ecrt.dat$mean), col="grey")
for (i in 1:10){
	lines(dem, probs.list[[i]][,2], col="black", lty=1, lwd=3)
	lines(supsmu(dem, probs.list[[i]][,1]), col="gray48", lty=1, lwd=3)
	lines(supsmu(dem, probs.list[[i]][,3]), col="gray48", lty=1, lwd=3)
	}
rug(ecrt.dat$mean[ecrt.dat$ecrt==1], lwd=2)
box()

plot(density(newest.dat$mean,na.rm=T), ylim=c(0, 1), xlim=c(-2, 2.5), main="ICCPR", xlab="Latent Democracy", ylab="Pr(Ratification)", axes=F, type="n")
abline(h=c(0,0.2,0.4,0.6,0.8,1), col=tr.gr, lty=3, lwd=1.5)
polygon(density(newest.dat$mean,na.rm=T), col="grey")
axis(1, at=c(-2, -1, 0, 1, 2))
axis(2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1))
polygon(density(newest.dat$mean,na.rm=T), col="grey")
for (i in 1:10){
	lines(dem, probs.list.2[[i]][,2], col="black", lty=1, lwd=3)
	lines(supsmu(dem, probs.list.2[[i]][,1]), col="gray48", lty=1, lwd=3)
	lines(supsmu(dem, probs.list.2[[i]][,3]), col="gray48", lty=1, lwd=3)
	}
rug(newest.dat$mean[newest.dat$iccpr_rat==1], lwd=2)
box()

dev.off()

#### effect of transition

probs.list<-list()
for (i in 1:50){

	probs<-matrix(nrow=30, ncol=4)

	for (j in 1:30){
		pr1<-
		plogis(
		ests.dif.list[[i]][,46]+
		ests.dif.list[[i]][,42]*0+
		ests.dif.list[[i]][,43]*0+
		ests.dif.list[[i]][,44]*j
		)
		
		pr2<-
		plogis(
		ests.dif.list[[i]][,46]+
		ests.dif.list[[i]][,42]*1+
		ests.dif.list[[i]][,43]*0+
		ests.dif.list[[i]][,44]*j
		)
		
		pr<-pr2-pr1
		probs[j,1]<-quantile(pr, probs=.05)
		probs[j,2]<-mean(pr)
		probs[j,3]<-quantile(pr, probs=.95)
		probs[j,4]<-j

		}

	probs.list[[i]]<-probs

	}

lab1<-expression(Democracy^2)
lab2<-expression(Delta[t-1])
lab3<-expression(Delta[t-2])
lab4<-expression(Delta[t-3])
lab5<-expression(Delta[t-4])
lab6<-expression(Delta[t-5])
dif.labs<-c(lab2,lab3,lab4,lab5,lab6)

probs.t25<-matrix(nrow=length(probs.list),ncol=4)
for (i in 1:length(probs.list)){
	probs.t25[i,]<-probs.list[[i]][20,]
	}

pdf("ecrt_dif_20.pdf")
par(mfrow=c(5,1), mar=c(1, 2, 2, 1), cex.axis=.8)
for (j in 1:5){
	plot(seq(1:10), probs.t25[((j*10)-(9)):(j*10),2], type="p", pch=16, cex=1.75, ylab="", xlab="" , 				ylim=c(-0.3,0.3), main=dif.labs[j], cex.main=1.5, axes=F)
	for (i in 0:9){
		lines(c(i+1, i+1), c(probs.t25[((j*10)-(9))+i,1],probs.t25[((j*10)-(9))+i,3]), lwd=1.75)
		}
	axis(1, at=seq(1:10), labels=c("","","","","","","","","",""))
	axis(2, at=c(-0.4, -0.2, 0, 0.2), labels=c("-0.4", "-0.2", "0", "0.2"))
	abline(h=0, lty=2)
	abline(h=c(-0.2,0.2), lty=3, col=tr.gr, lwd=1.5)
	box()
	}
dev.off()	

### Time series graphs
ecrt.dat$ccode[ecrt.dat$ccode==316&ecrt.dat$year<1993]<-315
ecrt.dat<-subset(ecrt.dat,ecrt.dat$ccode!=316)
cc<-c(230,235,350,290,315,355)
pdf("ecrt_ts1.pdf")
par(mfrow=c(3, 2), mar=c(3, 2, 4, 1))
for (i in 1:length(cc)) {
	plot(ecrt.dat$mean[ecrt.dat$ccode==cc[i]]~ecrt.dat$year[ecrt.dat$ccode==cc[i]], type="l", lwd=2, xlab="", ylab="Latent Democracy", main=ecrt.dat$country.y[ecrt.dat$ccode==cc[i]][1], ylim=c(-1.75, 2.5))	

	polygon(c(ecrt.dat$year[ecrt.dat$ccode==cc[i]],rev(ecrt.dat$year[ecrt.dat$ccode==cc[i]]),ecrt.dat$year[ecrt.dat$ccode==cc[i]][1]),c(ecrt.dat$pct975[ecrt.dat$ccode==cc[i]],rev(ecrt.dat$pct025[ecrt.dat$ccode==cc[i]]),ecrt.dat$pct025[ecrt.dat$ccode==cc[i]][1]),col=tr.gr,border=NA)

	abline(v=ecrt.dat$year[ecrt.dat$ccode==cc[i]&ecrt.dat$ecrt==1][1], lty=2)
	abline(h=c(-1,0,1,2), lty=3, col=tr.gr, lwd=1.5)
}
dev.off()

ecrt.dat$country.y<-as.character(ecrt.dat$country.y)
ecrt.dat$country.y[ecrt.dat$cc==200]<-"UK"
ecrt.dat$country.y[ecrt.dat$cc==260]<-"W. Germany"
ecrt.dat$country.y[ecrt.dat$cc==265]<-"E. Germany"
ecrt.dat$country.y[ecrt.dat$cc==365]<-"USSR/Russia"
ecrt.dat$country.y[ecrt.dat$cc==345]<-"Yugoslavia"
cc<-c(200, 205, 210, 211, 212, 220, 225, 260, 305, 325, 338, 375, 380, 385, 390, 395)
pdf("ecrt_ts2.pdf")
par(mfrow=c(4, 4), mar=c(3, 2, 3, 1))
for (i in 1:length(cc)) {
#	tmp<-paste(as.character(cc[i]),"pdf",sep=".")
#	pdf(tmp)
	plot(ecrt.dat$mean[ecrt.dat$ccode==cc[i]]~ecrt.dat$year[ecrt.dat$ccode==cc[i]], type="l", lwd=2, xlab="", ylab="Latent Democracy", main=ecrt.dat$country.y[ecrt.dat$ccode==cc[i]][1], ylim=c(-1.75, 2.5), axes=F)	

	polygon(c(ecrt.dat$year[ecrt.dat$ccode==cc[i]],rev(ecrt.dat$year[ecrt.dat$ccode==cc[i]]),ecrt.dat$year[ecrt.dat$ccode==cc[i]][1]),c(ecrt.dat$pct975[ecrt.dat$ccode==cc[i]],rev(ecrt.dat$pct025[ecrt.dat$ccode==cc[i]]),ecrt.dat$pct025[ecrt.dat$ccode==cc[i]][1]),col=tr.gr,border=NA)

	abline(v=ecrt.dat$year[ecrt.dat$ccode==cc[i]&ecrt.dat$ecrt==1][1], lty=2)
	abline(h=c(-1,0,1,2), lty=3, col=tr.gr, lwd=1.5)

	axis(1)
	axis(2, at=c(-1, 0, 1, 2), labels=c("-1", "", "", "2"))
	box()
#	dev.off()
}
dev.off()

cc<-c(265, 310, 339, 343, 344, 345, 346, 352, 359, 360, 365, 366, 367, 369, 370, 371, 372, 373, 640)
pdf("ecrt_ts3.pdf")
par(mfrow=c(5, 4), mar=c(3, 2, 3, 1))
for (i in 1:length(cc)) {
#	tmp<-paste(as.character(cc[i]),"pdf",sep=".")
#	pdf(tmp)
	plot(ecrt.dat$mean[ecrt.dat$ccode==cc[i]]~ecrt.dat$year[ecrt.dat$ccode==cc[i]], type="l", lwd=2, xlab="", ylab="Latent Democracy", main=ecrt.dat$country.y[ecrt.dat$ccode==cc[i]][1], ylim=c(-1.75, 2.5), axes=F)	

	polygon(c(ecrt.dat$year[ecrt.dat$ccode==cc[i]],rev(ecrt.dat$year[ecrt.dat$ccode==cc[i]]),ecrt.dat$year[ecrt.dat$ccode==cc[i]][1]),c(ecrt.dat$pct975[ecrt.dat$ccode==cc[i]],rev(ecrt.dat$pct025[ecrt.dat$ccode==cc[i]]),ecrt.dat$pct025[ecrt.dat$ccode==cc[i]][1]),col=tr.gr,border=NA)

	lines(ecrt.dat$pct975[ecrt.dat$ccode==cc[i]]~ecrt.dat$year[ecrt.dat$ccode==cc[i]], col="grey", lwd=2)
	lines(ecrt.dat$pct025[ecrt.dat$ccode==cc[i]]~ecrt.dat$year[ecrt.dat$ccode==cc[i]], col="grey", lwd=2)	

	abline(v=ecrt.dat$year[ecrt.dat$ccode==cc[i]&ecrt.dat$ecrt==1][1], lty=2)
	abline(h=c(-1,0,1,2), lty=3, col=tr.gr, lwd=1.5)

	axis(1)
	axis(2, at=c(-1, 0, 1, 2), labels=c("-1", "", "", "2"))
	box()
#	dev.off()
}
dev.off()


#### Ropeladder plots

### coefficient estimates from E Court models
dem.ests<-matrix(nrow=10,ncol=3)
demsq.ests<-matrix(nrow=10,ncol=3)
for (i in 1:length(ests.list)){
	dem.ests[i,1]<-mean(ests.list[[i]][,42])
	dem.ests[i,2:3]<-quantile(ests.list[[i]][,42],c(0.05,0.95))
	demsq.ests[i,1]<-mean(ests.list[[i]][,44])
	demsq.ests[i,2:3]<-quantile(ests.list[[i]][,44],c(0.05,0.95))
	}

dif.ests<-matrix(nrow=50,ncol=3)
for (i in 1:length(ests.dif.list)){
	dif.ests[i,1]<-mean(ests.dif.list[[i]][,42])
	dif.ests[i,2:3]<-quantile(ests.dif.list[[i]][,42],c(0.05,0.95))
	}

coefs<-matrix(nrow=9,ncol=3)
coefs[1,1]<-mean(dif.ests[41:50,1])
coefs[1,2]<-min(dif.ests[41:50,2])
coefs[1,3]<-max(dif.ests[41:50,3])

coefs[2,1]<-mean(dif.ests[31:40,1])
coefs[2,2]<-min(dif.ests[31:40,2])
coefs[2,3]<-max(dif.ests[31:40,3])

coefs[3,1]<-mean(dif.ests[21:30,1])
coefs[3,2]<-min(dif.ests[21:30,2])
coefs[3,3]<-max(dif.ests[21:30,3])

coefs[4,1]<-mean(dif.ests[11:20,1])
coefs[4,2]<-min(dif.ests[11:20,2])
coefs[4,3]<-max(dif.ests[11:20,3])

coefs[5,1]<-mean(dif.ests[1:10,1])
coefs[5,2]<-min(dif.ests[1:10,2])
coefs[5,3]<-max(dif.ests[1:10,3])

coefs[7,1]<-mean(demsq.ests[,1])
coefs[7,2]<-min(demsq.ests[,2])
coefs[7,3]<-max(demsq.ests[,3])
coefs[8,1]<-mean(dem.ests[,1])
coefs[8,2]<-min(dem.ests[,2])
coefs[8,3]<-max(dem.ests[,3])

### coefficient estimates from CCPR models 

dem.ests.2<-matrix(nrow=10,ncol=3)
demsq.ests.2<-matrix(nrow=10,ncol=3)
for (i in 1:length(ests.list.2)){
	dem.ests.2[i,1]<-mean(ests.list.2[[i]][,31])
	dem.ests.2[i,2:3]<-quantile(ests.list.2[[i]][,31],c(0.05,0.95))
	demsq.ests.2[i,1]<-mean(ests.list.2[[i]][,33])
	demsq.ests.2[i,2:3]<-quantile(ests.list.2[[i]][,33],c(0.05,0.95))
	}

dif.ests.2<-matrix(nrow=50,ncol=3)
for (i in 1:length(ests.dif.list.2)){
	dif.ests.2[i,1]<-mean(ests.dif.list.2[[i]][,31])
	dif.ests.2[i,2:3]<-quantile(ests.dif.list.2[[i]][,31],c(0.05,0.95))
	}

coefs.2<-matrix(nrow=9,ncol=3)
coefs.2[1,1]<-mean(dif.ests.2[41:50,1])
coefs.2[1,2]<-min(dif.ests.2[41:50,2])
coefs.2[1,3]<-max(dif.ests.2[41:50,3])

coefs.2[2,1]<-mean(dif.ests.2[31:40,1])
coefs.2[2,2]<-min(dif.ests.2[31:40,2])
coefs.2[2,3]<-max(dif.ests.2[31:40,3])

coefs.2[3,1]<-mean(dif.ests.2[21:30,1])
coefs.2[3,2]<-min(dif.ests.2[21:30,2])
coefs.2[3,3]<-max(dif.ests.2[21:30,3])

coefs.2[4,1]<-mean(dif.ests.2[11:20,1])
coefs.2[4,2]<-min(dif.ests.2[11:20,2])
coefs.2[4,3]<-max(dif.ests.2[11:20,3])

coefs.2[5,1]<-mean(dif.ests.2[1:10,1])
coefs.2[5,2]<-min(dif.ests.2[1:10,2])
coefs.2[5,3]<-max(dif.ests.2[1:10,3])

coefs.2[7,1]<-mean(demsq.ests.2[,1])
coefs.2[7,2]<-min(demsq.ests.2[,2])
coefs.2[7,3]<-max(demsq.ests.2[,3])
coefs.2[8,1]<-mean(dem.ests.2[,1])
coefs.2[8,2]<-min(dem.ests.2[,2])
coefs.2[8,3]<-max(dem.ests.2[,3])

ylabels1<-c(lab6, lab5, lab4, lab3, lab2)
lab0<-expression(Democracy)
ylabels2<-c(lab1, lab0)
ylabels4<-c("Models 2-6 ", " Model 1 ")

pdf("ecrt_coefs.pdf")
par(mfrow=c(1,2), oma=c(1, 4, 1, 1), mar=c(2, 1.5, 1.5, 1), cex.axis=0.85, cex.main=1.25)

plot(coefs[,1], seq(1:dim(coefs)[1]), type="p", pch=16, cex=1.5, ylab="", xlab="", axes=F, xlim=c(-4, 6), main="European Court")
for (i in 1:dim(coefs)[1]){
	lines(coefs[i,2:3], c(i, i), lwd=1.5)
}
axis(2, labels=ylabels1, at=seq(1:5), las=1)
axis(2, labels=ylabels2, at=c(7, 8), las=1)
axis(2, labels=ylabels4, at=c(6, 9), las=1, font=2, tick=F)
axis(1, at=seq(from=-4, to=6, by=2))
abline(v=0, lty=2)
abline(v=c(-4,-2,2,4,6), lty=3, col=tr.gr, lwd=1.5)
box()

plot(coefs.2[,1], seq(1:dim(coefs)[1]), type="p", pch=16, cex=1.5, ylab="", xlab="", axes=F, xlim=c(-4, 6), main="ICCPR")
for (i in 1:dim(coefs.2)[1]){
	lines(coefs.2[i,2:3], c(i, i), lwd=1.5)
}
axis(1, at=seq(from=-4, to=6, by=2))
abline(v=0, lty=2)
abline(v=c(-4,-2,2,4,6), lty=3, col=tr.gr, lwd=1.5)
box()
dev.off()
